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SECTION  I 
INTRODUCTION 

1.  GENERAL  REVIEW 

As  a  step  toward  the  realization  of  the  concept  of  a  "Numerical 

Wind  Tunnel",  the  unified  use  of  the  full  continuum  flow  equations 

(Navier-Stokes)  was  employed  to  yield  numerical  solutions  for 

supersonic  flow  past  nozzle-afterbody  configurations  for  arbitrary 

nozzle  geometries  and  flow  conditions.  This  particular  flow 

problem  is  usually  investigated  through  wind  tunnel  testing.  The 

complexity  of  the  flow  structure  and  the  lack  of  confidence  in 

existing  numerical  capabilities  are  among  the  factors  that  deterred 

investigators  from  attacking  the  problem  using  the  present  approach. 

1-3 

Numerical  approaches,  however,  were  tried  by  several  investigators 
using  patching  procedures  for  different  regions,  or  other  restrictive 
assumptions  such  as  solid  plume  simulators. 

This  work  was  done  for  the  purpose  of  developing  different 

techniques  for  estimating  the  aerodynamic  properties  of  aircraft 

components  as  a  step  forward  toward  the  design  integration  goal. 

This  work  is  the  second  stage  of  a  two-stage  project.  The 
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completed  work  of  the  air  intake  inlet  problem  was  the  first  stage. 

The  aerodynamics  of  propulsion  elements  in  aircraft  and 
rockets  plays  a  significant  role  in  determining  the  importance  of 
the  airborn  vehicle.  The  pressure  drag  of  the  jet  engine  exhaust 
nozzle,  for  example,  contributes  considerably  to  the  total  drag 
of  the  air  worthy  jet  propelled  vehicle.  Its  minimization,  there¬ 
fore,  receives  considerable  attention  by  experimentalists  for 
optimum  boattailing.  The  pressure  drag  on  the  nozzle  surface 
can  not  be  properly  computed  without  considering  the  mutual 
effects  of  the  exhaust  jet  and  the  outer  stream.  Therefore,  the 
domain  of  mutual  influence  between  the  outer  and  the  existing 
jet  flow  must  be  considered.  This  mutual  effect  is  being  con¬ 
sidered  by  the  present  approach. 
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This  particular  flow  problem  is  directly  related  to  the 
mixing  of  two  supersonic  axisymmetric  streams.  It  also  involves 
most  of  the  flow  features  that  can  be  observed  in  a  compressible 
supersonic  flow  field.  Turbulent  boundary  layer,  separation, 
free  reattachment  (not  on  a  surface) ,  lip  shock,  separation  shock, 
barrel  shock,  shock  intersection  and  reflection,  Mach  disk,  slip 
surface,  free  turbulent  mixing,  "inviscid"  plume  structure  and 
other  interesting  details  can  be  observed. 

The  basic  flow  structure  is  described  briefly  next. 

2.  FLOW  STRUCTURE 

Depending  on  the  nozzle  pressure  ratio  (NPR)  which  will  be 
defined  here  as  NPR  =  P^/P^  for  reasons  discussed  in  Section  II, 
the  flow  structure  can  assume  one  of  two  familiar  modes. 

Considering  only  underexpanded  nozzles,  and  if  the  NPR  is 
low  (lower  than  a  certain  critical  value  NPR  )  the  turning  angle 
of  the  exiting  flow  at  the  nozzle  wall  will  be  relatively  small, 
allowing  the  expansion  fan  to  reflect  at  the  center  line  r  =  0, 
with  incident  and  reflected  angles  almost  equal  to  90°.  This 
will  result' in  the  familiar  repeated  "diamond"  pattern,  as  shown 
in  Figure  1 . 

If  the  nozzle  boattail  outline  curve  is  of  steep  gradient, 
flow  separation  may  occur  on  the  surface,  however,  reattachment  is 
not  expected  to  occur  at  the  solid  wall,  but  rather  down  stream. 

The  separation  shock  may  coalesce  with  the  weak  attachment  shock 
forming  one  X  (lambda)  structure.  Two  circulating  flow  bubbles,  or 
more,  are  mostly  expected  depending  on  the  geometry  of  the  nozzle  lip 
and  the  NPR.  The  mixing  between  the  external  and  the  jet  streams 
will  take  place  along  the  "inviscid  plume  boundary"  which  disappears, 
giving  way  to  a  widening  mixing  "layer"  that  increases  in  "width" 
and  finally  intersects  with  the  line  of  symmetry. 

For  larger  NPR  (larger  than  the  critical  NPR)  the  jet  flow 
turning  angle  at  the  nozzle  lip  is  relatively  large,  allowing 
the  expansion  fan  to  reflect  at  the  center  line  with  angle 
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smaller  than  90°,  which  reflects  again  at  the  ficticious  "inviscid 
plume  boundary"  in  a  form  that  allows  them  to  coalesce  forming  one 
shock  known  as  the  barrel  shock  before  they  intersect  with  the 
center  line,  as  shown  in  Figure  (1.1B)  with  the  familiar  triple 
point  formation.  The  intersection  at  the  center  line  is  in  the 
form  of  a  curved  Mach  disk  with  subsonic  region  behind  it.  A 
slip  surface  initiates  at  the  triple  point  and  extends  behind  the 
Mach  disk. 

For  the  same  boattail  geometry,  flow  separation  on  the  boat- 
tail  surface  is  expected  to  start  earlier  upstream  for  the  high 
NPR.  In  addition,  the  large  flow  turning  angle  will  then  cause 
the  shock  formed  to  be  stronger.  The  viscous  interaction  is 
much  stronger  in  this  case  and  no  well-def ined  reattachment  point 
is  expected. 

3.  THE  PRESENT  APPROACH 

The  domain  of  present  interest  is  shown  in  Figure  (1) ,  where 
the  viscous  interaction  between  the  external  and  jet  streams  is 
a  dominant  factor.  Therefore,  the  use  of  the  complete  Navier- 
Stokes  equations  is  necessary  to  capture  such  interactions. 

The  equations  will  be  solved  uniformly  all  over  the  domain  of 
interest,  thus  avoiding  the  need  for  any  patching  procedure  as 
those  of  References  (3)  and  (5) .  This  uniformity  also  enables  the 
true  presentation  and  computation  of  the  mutual  effects  of  the 
jet  and  the  boattail  surface  flow  conditions  and  vice  versa,  with¬ 
out  the  need  for  superposition  of  effects  in  an  iterative  procedure 
to  correct  for  the  mutual  influences. 

A  coordinate  transformation  will  be  used  to  map  the  present 
domain  onto  a  square  of  unit  length.  The  mesh  points  in  that 
transformed  plane  are  uniformly  spaced,  while  the  corresponding 
points  in  the  physical  domain  are  highly  nonuniform  for  the  proper 
concentration  required  for  turbulent  flow. 

Since  the  present  case  of  interest  is  only  for  the  supersonic 

flow  at  M  =1.5,  the  numerical  scheme  chosen  should  be  a  shock 

00 

capturing  method  since  the  details  and  exact  location  of  ’the 
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shocks  is  not  known  a  priori.  Therefore,  MacCormack’s  time 
dependent,  explicit  scheme  will  be  favored  for  its  proven  reliability. 
Although  this  explicit  scheme,  as  most  explicit  schemes,  have 
severe  stability  restrictions,  it  is  to  be  understood  here,  that 
it  is  not  the  task  of  the  present  work  to  test  and  utilize  more 
efficient  schemes,  which  are  relatively  recent  leaving  that  to 
simpler  problems  for  future  investigation. 

The  present  work  is  considered  pioneering  due  to  the  very 
little  known  information  about  the  flow  conditions  in  the 
immediate  vicinity  of  the  computational  region.  In  addition, 
several  engineering  judgements  had  to  be  made  regarding  the  flow 
conditions,  shape  and  extension  of  the  region  of  computation,  the 
coordinate  system  chosen/  mesh  point  distribution  and  the  turbulence 
model  used. 

Therefore,  the  major  objective  in  this  work  is  to  show  that 
such  a  complex  flow  can  be  computed  successfully,  rather  than 
emphasizing  accuracy  or  computational  efficiency. 
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SECTION  II 

FORMULATION  OF  THE  PROBLEM 


In  this  section,  details  of  the  steps  of  formulation  are 
presented.  The  choice  of  the  test  case,  the  coordinate  mapping 
used,  the  form  of  the  governing  equations,  boundary  conditions 
specified,  turbulence,  model  and  other  pertinent  details  are 
discussed. 

1.  THE  SPECIFIC  TEST  CASE  CHOSEN 

To  test  and  verify  the  numerical  computation,  experimental 
results  must  be  available  for  comparison.  Limited  experimental 
results  are  available  for  the  almost  parallel  mixing  of  supersonic 
streams.  For  nonparallel  streams,  there  are  many  results.  For 
example,  results  for  single  jet  exhausting  into  static  air,  or  the 
mixing  of  exhaust  with  subsonic  external  stream. 

In  Reference  (6) ,  three  AGARD  convergent  nozzles  were  tested, 
but  only  the  boattail  surface  pressure  was  reported.  The  three 
nozzles  are  denoted  by  10°,  15°,  25°  nozzles  and"  were  chosen  by 
the  AGARD  Organization  to  study  the  deviation  in  the  measurements 
reported  by  the  different  testing  facilities  for  the  same  geometries. 
In  the  mentioned  reference,  measurements  for  cases  with  hot  jet 
exhaust  and  different  free  stream  Mach  number  (supersonic  and 
subsonic)  were  also  presented  for  some  of  these  three  geometries. 

No  other  flow  variables  or  information  down  the  nozzle  exit  or  on 
the  boattail  surface  were  provided. 

For  the  present  work,  the  10°-nozzle  was  chosen  with  Mach 
number  MMe  =  1.5,  at  NPR  =  7.09.  The  exact  geometry  of  that 
nozzle  is  given  in  Figure  (2.1).  The  experiments  on  that  nozzle 
were  run  in  the  AEDC  supersonic  wind  tunnel,  with  Reynolds  number 
=  2.5  x  10  per  foot.  .  At  the  station  130.47"  from  the  tip  nose 
of  the  engine  model,  the  local  Reynolds  number  is  then  27.18  x  10p 
signifying  fully  turbulent  flow  conditions. 
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Figure  2.  Geometry  of  the  AGARD  10°-Nozzle. 
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Station 

157.121 


Figure  3.  The  AEDC  Nozzle  Boattail  Geometry. 
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Another  nozzle  geometry  was  used  by  the  AEDC  at  the  early  stage 
of  its  study,  however  no  experimental  data  is  available  for  super¬ 
sonic  external  flow.  This  nozzle  will  be  denoted  here  as  the  AEDC 
early  model.  This  geometry  was  used  in  the  present  work  to 
illustrate  the  effect  of  the  boattail  geometry  on  the  surface  c  - 
distribution  in  comparison  with  the  result  for  the  AGARD  10° 
nozzle.  Details  of  this  nozzle  geometry  are  given  by  Figure  (2.2). 

It  is  to  be  noted  that  the  value  7.09  for  the  NPR  as  defined 

here  by  P j/Pooe  is  not  "very"  high.  Usually  the  NPR  is  defined  in 

the  case  of  a  right  cylinder  with  no  external  flow,  so  as  to  reflect 

the  ratio  between  the  static  pressure  value  at  the  exit  plane,  P j  , 

and  the  ambient  value  in  the  immediate  vicinity  of  the  nozzle  lip, 

P  ,  as  demonstratedv in  Figure  (2.3  a,b) .  However,  for  curved 
a 

boattail  walls  with  external  flow,  the  value  of  the  static  pressure 

at  the  lip,  P  ,  as  of  Figure  (2.3a)  is  not  known  a  priori, 
a 

Besides,  P  at  the  nozzle  lip  is  greatly  different  and  smaller 

(for  supersonic  flow)  than  the  value  of  the  undisturbed  pressure 

of  the  external  flow,  P  . 

e 

Therefore,  in  this  case,  since 

the  NPR  cannot  be  defined  as 

Pj/P&f  it  is  usually  defined  as 

P./Pn  .  For  M  =1.5,  the 
j '  °e  °°e 

ratio  P./P  is  only  3.7  for 
j  °°e  x 

the  corresponding  ratio  P^/ 

Poe  of  value  7.09.  Finally, 

it  is  to  be  noted  that  the 

critical  value  NPR  =  2.  - 

cr 

2.4,  is  only  valid  for  the 

case  (b)  of  Figure  (2.3), 

where  here  NPR  denotes  P ./P  . 

T  a 


Figure  4.  Nozzle  Pressure  Ratio 
Definition. 
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2. 


COORDINATE  SYSTEM 


Domain  of  Computation 

The  cylindrical  coordinates  (r,  0,  x)  were  used  for  this 
axisymmetric  geometry  where  r  is  in  the  radial  direction  and  x  is 
along  the  axis  of  symmetry.  The  0  dependency  is  dropped  for  the 
axisymmetry  formulation.  Of  the  16.0  inches  of  the  nozzle  body 
shown  in  Figure  (2.1),  only  10.0  inches  were  considered  to  minimize 
the  number  of  points  needed  for  computation.  A  height  of  12.0 
inches  was  taken  normal  to  the  boattail  surface  representing  about 
four  times  the  estimated  boundary  layer  thickness  at  that  location. 
Along  the  axis,  16.0  inches  representing  about  four  times  the  exit 
diameter  was  considered  appropriate,  ensuring  that  the  "far  field" 
boundary  conditions  will  be  in  the  supersonic  region  if  a  Mach 
disc  followed  by  the  subsonic  region  should  occur  at  higher  NPR. 

Coordinates  Transformation 

< 

The  (r,  x)  physical  domain  of  Figure  (2.4)  was  mapped  to  a 
unit  square  in  the  computational  plane  (ri,£)  through  the  mapping 
procedure  described  in  Reference  (7) .  The  n  =  Constant  lines  are 
aligned  parallel  to  the  boattail  wall  surface  and  the  £  =  constant 
lines  are  in  the  direction  normal  to  the  boattail  surface.  Along 
the  n  =  constant  lines,  39  points  were  used,  and  30  points  were 
utilized  along  the  £  =  constant  lines. 

It  was  found  difficult,  as  will  be  discussed  next,  to  obtain 
the  proper  source  distribution  necessary  for  the  procedure  of 
Reference  (7)  that  would  cause  the  coordinate  lines  to  satisfy 
all  the  requirements.  Therefore,  an  averaging  procedure  was  used, 
with  a  sacrifice  in  the  total  smoothness  of  the  transformation 
coefficients  r\  ,  n  /  £  and  £  (subscripts  denote  differentiation 

27  X  3T  X 

with  respect  to  the  subscript  variable) . 

In  the  mentioned  mapping  procedure,  the  point  distribution 
should  be  provided  along  the  boundary  afhjk  of  Figure  (2.4).  The 
appropriate  points  along  af  as  well  as  the  corresponding  points 
along  hj  were  chosen.  Along  the  boundaries  fh  and  a j ,  the 
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distributions  were  basically  similar,  both  with  exponential  growth 
allowing  the  first  point  to  be  about  0 . 07  inches  from  the  wall  and 
the  last  point  12.0  inches  away.  Although  0.07  inches  is  considered 
very  large  with  regard  to  turbulent  flow,  it  was  used  here  to  allow 
relatively  large  time  step  for  the  time  dependent  method  as  will 
be  discussed  later  in  this  section.  This  fact  will  therefore  lead 
to  large  inaccuracies  in  the  estimated  values  for  the  viscous 
drag  coefficient  cf.  However  c^,  which  is  the  main  parameter 
of  interest  here  ,  is  not  greatly  sensitive  to  the  mesh  step  size 
as  long  as  it  is  small  enough. 

Requirements  in  the  Chosen  Coordinates 

The  following  requirements  were  set  for  determining  the 
appropriate  coordinate  system.  They  also  represent  a  difficulty 
in  determining  the  proper  source  distribution  for  the  method  of 
Reference  (7) . 

1.  Maintaining  orthogonality  to  the  boattail  surface,  especially 
near  the  surface,  to  simplify  the  application  of  the  derivative 
boundary  conditions.  This  had  to  be  relaxed  in  the  region  of 

j 

Figure  (2.4)  to  allow  smooth  distribution  of  the  corresponding 
points  on  the  boundary  i j . 

2.  Orthogonality  is  also  desirable  near  the  outer  boundary 
hj  for  simple  application  of  the  "no  change"  boundary  condition. 

This  requirement  was  relaxed  locally. 

3.  More  dense  points  are  needed  in  the  region  be,  compared 

to  ab,  so  that  flow  separation  can  be  detected  and  numerical 

instability  is  possibly  avoided. 

4.  In  region  cf,  a  small  step  size  is  needed  near  the  nozzle 
lip  in  the  region  cd,  however,  larger  steps  can  be  allowed  in 
region  de.  Relatively  small  steps  in  the  region  ef  are  needed 
near  the  center  line  to  minimize  the  error  of  the  described 
derivative  boundary  conditions  especially  when  only  first  order 
accurate  differences  are  used  at  the  boundaries. 
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5.  More  points  are  needed  near  the  wall,  in  the  region  ak 
to  capture  the  large  gradient  of  the  turbulent  boundary  layer. 
However,  if  the  coordinate  lines  n  -  constant  are  parallel  to  the 
wall  af,  it  will  be  very  difficult  to  get  this  distribution  to 

be  different  in  the  region  fg  where  not  as  many  points  might  be 
required.  In  fact,  even  for  jet  flow  exiting  at  =  1.001,  with 
lines  of  characteristics  almost  normal  to  the  center  line  fh,  no 
difficulty  was  encountered  with  relatively  wide  spread  points. 

6.  The  coordinate  lines  in  the  region  between  df  and  ih 
should  be  as  orthogonal  as  possible  to  the  center  line  fh,  to 
simplify  the  implementation  of  the  boundary  conditions. 

Remarks  and  Comments  on  the  Present  Coordinates 

The  sharp  corner  "A"  of 
nozzle  lip  was  first  avoided 
by  placing  the  points  (I  =  24, 

25)  on  each  side  as  shown  in 
Figure  (2.5A).  This  was  done 
to  have  continuous  trans¬ 
formation  derivatives.  How¬ 
ever,  it  was  found  essential, 
later,  to  have  the  jet  exit 
plane  and  the  lip  line  be 
substituted  by  a  curve  AB  as 
that  of  Figure  (2.5B)  .  This 
curve  was  chosen  arbitrarily  and 
for  the  numerical  procedure. 

Although  placing  the  boundary  conditions  on  that  curve 
instead  of  the  original  straight  line,  will  not  represent  a  large 
approximation  as  will  be  seen  in  Section  (2.5),  it  is  undoubtedly 
an  undesirable  parameter.  The  possible  interpretation  will  also 


represents  an  additional  parameter 


be  discussed  in  Section 
(2.5).  The  exact  geometry 
of  the  back  shoulder  is 
presented  in  Figure  (2.6). 

As  was  mentioned 
earlier,  two  different 
nozzle  geometries  were 
examined.  The  AEDC  early 
model  was  first  utilized 
and  the  grid  point  distri¬ 
bution  along  the  wall  is 
shown  in  Figure  (2.6). 

When  the  AGARD  10° -nozzle 
was  used,  a  different  and 
more  uniform  point  distri¬ 
bution,  as  shown  by  curve  A 
in  Figure  (2.7)  was  attempted.  However,  this  new  distribution 
together  with  the  resulting  transformation  coefficients  led  to  blow 
ups  during  computations  due  to  negative  physical .quantities  near 
the  nozzle  lip.  Surprisingly,  this  problem  was  alleviated  when 
the  point  distribution  for  the  new  geometry  was  chosen  to  be  affine 
to  the  distribution  for  the  AEDC  model.  By  affine  it  is  meant  that 
only  the  r-coordinate  changes,  keeping  the  x-location  of  the  points 
the  same.  Both  the  affine  and  the  more  uniform  (but  not  affine) 
distributions  are  shown  in  Figure  (2.7).  The  possible  interpretation 
for  this  undesirable  behavior  is  that  the  affine  distribution 
have  more  points  near  the  nozzle  lip,  where  large  gradients  in 
pressure  and  density  occur. 

3.  THE  GOVERNING  EQUATIONS 

The  Navier-Stokes  equations  for  axisymmetric  and  two-dimensional 
turbulent  flows  can  be  written  in  the  following  familiar  form, 
where  the  dependent  variables  u,  v,  e  are  mass-averaged  as  described 
in  Reference  (8) ,  while  p  and  P  are  the  mean  (time  averaged) 
state  variables. 


Figure  6b.  Mesh  Points  Near  the 
Nozzle  Lip. 
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Figure  7.  Placement  of  Nozzle  Exit  Conditions. 
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Figure  8.  The  Affine  Distribution. 
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(2.4) 
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The  perfect  gas  relation  was  used  for  air,  i.e. 

P  =  pRT  (2.5) 

and  the  coefficient  of  viscosity  was  assumed  to  vary  with  tempera¬ 
ture  according  to  Sutherland's  law 

V  =  (2.27  X  10-8)  -(T  +T19'8'.'6)  (lbf  "  sec/ft2>  (2.6) 


The  laminar  and  turbulent  Prandtl  numbers  are  assumed  constant 
and  of  values  Pr  =  0.72  and  P  =  0.9  respectively. 

The  second  coefficient  of  viscosity  was  chosen  as 

At  =  (p  +  e)  •  6  (2.7) 


where  3  is  a  constant  less  or  equal  to  +1.0  (usually'  it  is 
negative  and  of  order- 10 ).  It  can,  of  course,  vary  as  a  function 
of  (r,  x)  .  It  is  used  during  the  early  stage  of  computation  to 
help  damping  large  transient  pressure  gradients  as  described  and 
used  in  Reference  (9)  .  However,  after  the  numerical  solution 
stabilizes  itself,  3  should  be  set  back  to  +1.0  to  represent  the 
true  Navier-Stokes  equations.  3  is  referred  to  as  the  normal 
stress  damping  coefficient  because  although  it  influences  particle 
friction,  it  only  occurs  in  the  normal  stress  components.  In 
addition,  it  is  only  effective  when  V  •  U  is  relatively  large, 
which  is  the  case  across  shocks  and  large  pressure  gradients  in 
general . 


For  numerical  computation,  Equation  (2.1)  is  written  in  the 
transformed  plane  of  (n-£)  as 


+  r 

at  L 


5  +  4-  f 

35  rj  + 

r  o 


a(r 


'joG)  ] 

^  J 


I  OF  ,  1 

+  K  a^+-T 


rJo 


(rjoG)  1  =  . 

an  J  Jo  • 


H 


r3o 


(2.8) 


19 


where  now  n  and  g  are  the  independent  variables  and  £  ,  £  ,  n  , 

x  r  x 

and  nr  are  the  four  transformation  derivatives  obtained  numerically 
from  the  mapping  procedure.  In  addition,  it  is  expected  that 
solving  these  equations  numerically  is  more  difficult  than  that  of 
Equation  (2.1)f  due  to  the  stiffness  introduced  by  the  stiff  trans¬ 
formation  derivatives  for  the  present  geometry. 


Equations  (2.8)  are  in  weak  conservative  form  due  to  the 
source  term  arising  only  for  the  axisymmetry  formulation,  and  due 
to  the  varying  coefficients  in  front  of  the  derivatives.  These 
equations  can  be  reduced  to  a  better  conservation  form  for  the 
axisymmetry  case  and  to  the  strong  conservation  form  for  the 
corresponding  2-D  case,  if  the  dependent  variables  as  regrouped 
according  to  Viviand"^ ,  the  form  of  the  equation  would  then  be 


3 

8t 


(2.9) 


where  J  is  the  Jacobian  of  the  transformation  defined  as 
c 

Jc  =  ?xnr  “  ?rV  , 

Although  this  form  was  not  used  in  the  present  work  it  is 
recommended  for  future  formulations  for  better  shock  capturing 
ability  and  possibly  better  accuracy. 


4.  TURBULENCE  MODELS 

The  experimental  tests  for  the  nozzle  geometry  under  considera¬ 
tion  (the  AGARD  10°-nozzle)  were  run  at  Reynolds  number  of 
2.5  x  10  /ft.  Therefore,  the  flow  is  expected  to  be  fully  turbulent 
with  a  Reynolds  number  in  excess  of  25.0  x  10  at  the  end  of  the 
long  cylindrical  forebody  at  station  130.47  inches  downstream  of 
the  nose  tip.  This  station  marks  the  beginning  of  the  nozzle 
boattail,  6.49  inches  upstream  of  the  station  where  computations 
begin  (see  Figure  (2.2)). 
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Complicated  eddy  viscosity  models  can  be  found  in  the 
literature  for  many  flow  patterns  including  the  axisymmetric  wake 
problem.  It  has  not  yet  been  demonstrated  that  the  use  of  sophisticated 
or  high  order  closure  models  would  necessarily  yield  better  results  in 
general.  Many  models  for  the  free  shear  layers,  jets  and  wakes 
using  the  turbulent  kinetic  energy  and  other  concepts  can  be 
found  for  example  in  References  (11)  and  (12) . 

It  was  decided  here  to  use  the  simplest  model  in  the  complex 
mixing  region.  The  eddy  viscosity  model  in  the  computational 
domain  was  split  into  different  models  in  different  regions  as 
shall  be  explained. 

The  computation  domain  is  split  into  three  regions  as  shown 
in  insert  A  of  Figure  (3).  In  region  I,  the  simple  two-layer 
algebraic  model  of  Reference  (13)  for  flat  plate  is  used  here 
since  the  effect  of  axisymmetry  was  estimated  to  be  minor  at  this 
flow  condition  for  a  cylinder  with  a  diameter  of  approximately 
4.0  inches.  The  inner  region  for  that  model  is  formulated  as 


ei  = p 
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and  the  outer  region  model  is  formulated  as: 
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where : 


r  is  the  normal  distance  measured  from  the  boattail  wall 
n 

x  is  the  wall  shear  stress 
w 

U.  is  the  tangential  velocity  in  direction  parallel  to  the 
wall 

U  is  the  average  tangential  velocity  outside  the  boundary 

layer  of  the  velocity  profile  at  that  station,  here 

taken  as  U,  at  the  point  of  index  j  =  15  (approximately 

r  *  2.0") 
n 


* 
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In  region  III,  the  eddy  viscosity  for  a  jet  is  used  which 
increases  in  the  axial  direction  while  remaining  constant  in  planes 
normal  to  the  axis.  The  following  model  was  used: 


'III 


-  e.(l  +  S) 

m 


where 

S  is  the  dimensionless  distance  along  the  "rays"  of  the 

coordinates  g  =  constant,  measured  from  the  jet  exiting 
station.  (S=0  at  the  exit  plane  and  S=1  at  the  far- 
field  boundary) 

e^n  is  the  initial  value  at  the  jet  exiting  station  (taken 
1  here  as  100  times  the  value  of  the  linear  viscosity 
coefficient  at  that  location) . 


In  region  II,  linear  variation  (with  respect  to  the  indexing 
parameter  i)  for  the  values  between  the  last  "ray"  in  region  I  and 
the  first  ray  of  region  III.  This  was  used  to  avoid  sudden  change 
in  the  e-value,  and  also  to  cover  the  region  of  possible  flow 
separation,  where  the  proper  turbulence  model  is  still  needed  to 
be  investigated. 

The  second  composite  model  of  sketch  B  of  Figure  (2)  was  also 
used,  but  its  results  are  not  presented  here.  No  major  differences 
were  observed  in  the  results  using  these  two  models.  In  this 
case,  region  III  of  model  A  is  split  further  into  two  separate 
regions,  III  and  IV.  In  the  new  region  III,  a  model  similar  to 
that  of  Reference. (14)  is  used  with  minor  alteration.  In  the 
present  work  the  model  used  was 


_1_  b  | 

eIII  RT  a  |Ui,  j  uq_,  J 


where 


b  is  the  "width"  of  the  mixing  region  as  shown  in 

Figure  (2) 

RT  is  the  Turbulent  Reynolds  Number  taken  here  as 

(390.  -  333.  e"(  -495  Mi, j) ) 

u  ..,T  ,  is  the  u-velocity  at  the  center  line  at  the  point  with 
U  '  the  index  j 


a 


is  a  factor  for  axisymmetry  effects 
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An  a  value  of  3.0  was  found  appropriate  for  the  present  case. 
Finally,  the  value  of  e  in  region  XV  was  set  to  zero,  representing 
the  inviscid  case  region.  The  value  of  e  in  region  II  was  also 
taken  as  the  linear  interpolation  between  the  values  of  regions 
I  and  III  as  described  previously. 

5.  BOUNDARY  CONDITIONS 

The  conditions  on  the  boundaries  of  the  physical  domain,  and 
the  corresponding  computational  domain  are  shown  in  Figure  (3) , 
and  they  are  discussed  next. 

The  Incoming  Flow  Conditions  at  Station  x  =  0  (g  =  0) 

The  Mach  number  and  all  the  flow  variables  are  not  known  at 
this  station,  therefore  some  engineering  judgement  was  required 
bearing  in  mind  the  major  objectives  of  the  present  work,  as 
stated  earlier.  It  is  assumed  that  the  Mach  number  is  still  1.5 
at  the  end  of  the  cylindrical  forebody  at  station  130.47  inches,, 
i.e.  that  the  compression  effect  at  the  nose  tip  is  exactly  offset 
by  the  expansion  at  the  forebody  shoulder  (see  JReference  (6)  for 
the  details  of  the  test  model  configuration) .  Moreover,  since 
computations  for  the  curved  boattail  start  at  station  6.49  inches 
down  from  station  130.47  inches,  the  use  of  the  value  M  =  1.5  is 
not  representative  of  the  real  flow.  In  fact,  the  tangential  Mach 
number  at  the  station  x  =  0  was  taken  as  1.74  based  on  a  Prandtl- 
Meyer  expansion  for  a  7°  turning  angle  for  the  corresponding  2-D 
flow. 

The  U  and  V  velocity  components  were  then  computed  at  this 
station  knowing  the  slope  of  the  body  at  this  location.  The 
turbulent  profile  for  the  tangential  velocity  was  assumed  to  be 
of  the  (1/7)  power  profile,  with  the  boundary  layer  thickness 
estimated  as  1.12  inches. 

The  static  temperature  and  pressure  profiles  were  assumed 
uniform  along  that  boundary,  with  values  of  0.905  and  0.7  Proe 

respectively  corresponding  to  the  7°  expansion  turn. 
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The  density,  the  specific  total  internal  energy  e,  and  the 
eddy  viscosity  profiles  were  then  computed  (using  the  equation  of 
state,  the  definition  of  e,  and  the  two  layer  turbulence  model)  and 
imposed  as  the  boundary  values. 

Conditions  on  the  Symmetry  Line  r  =  0  (£  =  1.0) 

Due  to  the  symmetry,  the  derivatives  in  the  direction  normal 
to  the  center  line  was  taken  to  be  zero  for  certain  variables. 

Since  not  all  the  variables  should  be  extrapolated  to  the  axis 
of  symmetry  in  this  manner,  otherwise  inconsistency  and  overspecifi¬ 
cation  would  result,  the  following  sequence  was  followed: 


.  3(pe)  3p  3u  3e  ,  ,  .  .  ,  . 

a)  — 5-= — —  ,  ttf,  ttf,  v  and  (pv)  were  set  to  zero,  then 

dt,  d  c,  d  c,  d  c, 

b)  P  was  computed  from  the  equation  of  state  using  p  and  e,  then 

c)  pu  was  computed  using  p  and  u. 

Conditions  on  the  Far-Field  Boundary  (n  =  1.0) 


The  variation  along  the  £  =  constant  lines  near  ri  =  1 . 0  was 
assumed  to  vanish.  Therefore,  the  following  conditions  were 
applied: 

v  3  ( pe )  3 p  3u  3e  3v  .  . 

a)  ITT'  3n'  an'  an  and  3n  were  set  to  zero'  then 

b)  P  was  computed  from  the  equation  of  state  using  p  and  e,  then 

c)  pu  and  pv  were  computed  using  p,  u  and  v. 

Conditions  at  the  Wall  and  at  the  Nozzle  Exit  (n  =  0. 0) 


On  the  boattail  wall: 


a)  u  and  v  are  zero  for  no  slip,  then  (pu) ,  (pv)  and  also  e 

are  zero.  T  is  set  at  T  ,  taken  equal  to  the  stagnation  temperature 

of  the  external  free  stream  TQe ,  to  simulate  the  wind  tunnel  tests 

or  T  can  be  specified  independently  for  hot  wall  cases 
w 

3P 

b)  ^  is  set  to  zero,  then 

c)  p  is  computed  from  the  equation  of  state  utilizing  P  and 
T ,  then 

d)  (pe)  is  computed  using  p,  u,  v  and  T 
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For  the  Nozzle  Exit  Plane 

a)  Uniform  profiles  u  =  u ^ ,  T  =  T  ^  ,  P  -  P..  and  v  ~  0  were 
imposed,  where  u ^ ,  1\  ,  P^  are  computed  from  the  given  experimental 
conditions  for  T0j,  poj'  Bej  an^  assuming  exit  Mach  number  of  1.01. 
Rej  is  the  jet  Reynolds  number  which  was  varied  in  the  experiments 
to  yield  the  required  NPR. 

b)  p  is  computed  from  the  equation  of  state,  and  e  at  this 
boundary  was  estimated  as  (100  p ^ ) 

c)  (pu) ,  (pv)  and  (pe)  were  then  computed  using  p,  u,  v  and  T. 

The  u  velocity  at  the  first  grid  point  in  the  jet  flow  nearest 
to  the  inner  nozzle  surface  was  set  to  0.5  u^  to  reflect  the  effect 
of  viscosity  and  simulate  the  corresponding  very  thin  boundary 

i  15 

layer 

As  mentioned  previously, 
a  smooth  curve  was  used  to 
join  the  points  AB  in  Figure 
(2.8),  instead  of  the 
vertical  line  AC.  This  will 
not  be  of  any  major  conse¬ 
quence  regarding  the  correct 
location  for  placing  the  cor¬ 
rect  boundary  conditions,  if  ment  for  the  Jet  Flow, 

the  line  AB  is  a  characteristic  line  for  the  supersonic  jet.  Then, 
according  to  the  linear  inviscid  supersonic  theory,  the  flow 
conditions  at  line  AC  remain  the  same  and  equal  to  those  at  line 
AB.  With  the  nozzle  assumed  in  chocking  condition,  M^  was  taken 
1.01  with  a  corresponding  characteristic  angle  a  =  81.9°. 

However,  line  AB  was  taken  at  another  angle  a  =  74°  for  smooth 
curvature.  This  difference  in  location  is  not  expected  to  affect 
the  solution  considerably,  especially  if  one  considers  the 
objectives  of  the  present  work. 

Although  the  smooth  shoulder  line  was  essential  for  obtaining 
successful  solutions,  its  geometry  is  not  expected  to  influence 
the 

gate  this  behavior. 
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SECTION  III 

THE  NUMERICAL  PROCEDURE 


The  time-dependent,  predictor-corrector  MacCormack 1 s  explicit 
16—18 

scheme  is  used  to  solve  Equation  (2.8)  as  successive  enhance¬ 

ment  in  time  toward  the  steady  state  flow  solution. 


1.  SPLITTING  PROCEDURE 


The  solution  is  obtained  by  two-step  splitting  in  the  r\-  and 
£-  directions.  Solution  advancement  in  time  can  start  in  the  n- 
and  then  the  £-  direction  or  vice  versa,  depending  on  which  is  larger, 
At^  or  At^  (as  defined  and  estimated  in  Section  (3.2)). 


of 


The  operator  L^(At^)  denotes  the  enhancement  of  the  solution 


3U  ,  3F  3  (r]°G) 

3t  4  3$  jQ  35 


0 


(3.1a) 


in  the  ^-direction  by  time  increment  of  At^  seconds.  Similarly, 

L  (At  )  represents  the"  similar  meaning  for  the  enhancement  of  the 
solution  of 


3U 

at 


3F  '  1 

nx  W  +  T 


Mr  °G) 

3n 


in  the  n-direction  by  At  increment. 


(3.1a) 


The  dependent  vector  variable  U(n,  t)  was  then  advanced 
in  time  as : 

D(n,?,At)  =  [l“/2  (At?/2)L  (MAt?)  L^/2(At?/2)]  •  U(n,?,t) 

with  At  =  MAt r  if  .Atw  <  At 

K  Z  n 


or  as 


U(n/5/At) 
with  At  =  NAt 

n 


[L 


”/2(Atn/2)  L^NAt^)  L 
if  At  <  At- 

n  5 


N/2 

n 


(Atn/2)] 


u(n,5,t) 

(3 . 2a ,b) 
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where  M  and  N  are  the  smallest  even  integers  of  the  quotients 
At  /Atr  and  Atr/At  ,  respectively  and  At  ,  Atr  are  the  maximum 
allowable  time  steps  in  the  n  and  £  directions  as  determined  by 
the  Courant-Friedrichs-Lewy  (CFL)  limit  with  viscous  stability 
requirement  as  will  be  shown  and  given  by  Equation  (3.3). 

The  values  of  M  and  N  were  chosen  to  be  close  to  2  for  the 
present  work  by  changing  the  grid  distribution.  This  fact  allows 
a  truly  alternating  direction  procedure  and  is  favored  for  the 
present  set  of  coordinates  to  enforce  the  boundary  conditions 
(especially  at  the  jet  exit  plane)  through  frequent  alterna¬ 
tion.  The  present  case  is  dissimilar  to  the  simple  flat 
plate  case  where  one  direction  (the  streamwise)  has  very  little 
gradients  and  the  need  for  frequent  alternation  of  the  solution 
path  is  unnecessary. 

The  maximum  time  step  as  determined  by  the  CFL  condition 
is  not  necessarily  very  near  the  wall,  but  rather  in  the  viscous 
mixing  region.  In  addition,  the  ratio  (At^/At^)  can  change 
drastically  from  less  to  greater  than  1.0  depending  upon  the 
relative  flow  conditions  of  the  external  flow  to"  the  jet  flow,  as 
was  noticed  for  the  different  flow  case  considered.  Therefore, 
form  3.2a  or  b  can  then  be  selected  accordingly.  This  procedure 
was  automated  in  the  computer  program. 


2.  STABILITY  CONDITION,  ESTIMATED  MAXIMUM  TIME  STEP 

-s. 

The  numerical  stability  requirement  for  MacCormack's  explicit 

scheme  is  governed  by  the  (CFL)  conditions  expressed  in  the  case 

17 

of  Cartesian  coordinates  as: 


At  =  minimum  (At  ;  At  ) 
•  *  x  y 


=  minimum 

i*  j 


Ax 

ui ,  j 

+  ci,j  +  max  if  <f 

r  rt 

&y 

1 

Ax' 

V-Vy+e)  l 

Ay  J 

|vi,j 

|  +  o  +  max  {  Y  (pP  +  p  ) 

^  -V*  V* 

1 

Ay' 

Ax  r 
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In  the  present  case, 
an  estimate  for  the 
maximum  time  step  was 
calculated  as: 

At  =  minimum  (At  ;AtJ 
if  j  n  * 


where 


Figure  10.  Nomenclature  for 

Estimating  the  Maximum 
Time  Step. 


where 


as?  =  y 

(xi,j  " 

i2  + 

(ri,j 

~  ri-lf jJ 

2 

1  / 

~ 

ASn  =  1 

,(xi,j  "  xi,j-i] 

77 

(r .  . 
if3 

1  / 

ur  =  Max .  (u.  .,v.  .,u.  ,  ,,v.  ,  .) 

5  ifD  if]  i-lfD  i-l rD 


Max.  (u.fj,  v  ,  u±fj_lf  v.^^) 


The  actual  time  step  finally  used  was  less  than  this  estimated 


maximum.  A  factor  of  0.7  was  used. 


3. 


NUMERICAL  DAMPING 


The  fourth  order  pressure  gradient  damping  concept  as  introduced 
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by  MacCormack  and  Baldwin  is  applied  here,  though  in  slightly 
different  form  to  suit  the  present  case.  The  constant  aQ  =  0.5 
for  the  flat  plate  case  as  used  in  that  reference,  for  example, 
has  now  been  made  a  function  of  the  location  (space) .  It  was  also 
logical  to  correlate  this  function  to  the  coordinate  transformation 
derivatives  n  ,  ri  r  and  £  . 

JL  2s.  XT 

This  damping  was  applied  in  both  computational  sweeps  in  the 
n-  and  £-  directions.  Some  numerical  experimentation  led  to  the 
choice  of  the  eight  constants  introduced  in  front  of  the  trans¬ 
formation  coefficients  to  allow  local  control  on  the  damping  in 
the  present  coordinates. 

This  damping  is  implemented  by  replacing 


F .  .  .  by  F .  .  .  +  F_ 

ll, D  ll,D  D±ifj 


and 


G.  .  .  by  G .  .  .  +  G_. 

11,3  11,3  D..fj 


in  both  the  predictor  and  corrector  steps  of  the  ^-solution  sweep 
of  Equation  (3.1a).  Here,  ii  =  i  and  i  +  1  for  the  predictor  and 
the  corrector  steps,  respectively. 


The  following  expressions  for  FQ  and  GD  were  used 


F 


=  cri  (  !  u .  .  .  I  +  c  .  .  . )  , 


Dii,j  ei  '  “'i 


P.  ,,  .  -  2P.  .  +  P.  ,  .1 

,  1+1,1 _ 1,3  1-1,3' 

ii, j  '  (P.^,  .  +  2P.  .  +  P.  T  .) 

,J  i+l,3  i,1  1-1,1 


3n  =ci  (  |  v .  .  .|+c..  .), 

Dii,j  nl 


P. .  . ,  .  -  2P.  .  .  +  P.  .  .  . 

,  n,1+l _ n,3  11,3-1' 

ii, j 7 '  (P. .  . +  2P. .  .  +  P.  .  .  , ) 

,J  n,3+l  ii,1  n,3-l 


where 


C_  =  0.02  ?2  +  0.12  ?2 

51  Xii,j  rii,3 


c  =  o.oi  n2  +  0.09  n2 
n,  r .  .  .  x .  .  . 

1  Ilf]  Ilf] 
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Then,  in  the  n-*sweep,  the  same  procedure  is  similarly  implemented 
by  replacing 


and 


F .  .  .  by  F .  .  .  +  F_ 

i ,33  y  i,33  D.  .. 

±  r  J  J 


G .  .  .  by  G .  .  .  +  G 

i/DD  iflD  D.  .. 

1  *  J  J 


where 


P. . ,  . .  -  2P.  .  .  +  P. 


and 


?  =  C  ( I u  .  I  +  c  )  1+1 'JJ _ Ll22 _ i-l,33 

D.  ..  '  1  i.r  j  j  1  i,jj  (P.  .n  .  .  +  2P.  .  .  +  P.  :  .  .) 

1,33  2  ,JJ  ,JJ  i+l  ,33  1,33  i-l ,33 


Ip  -  2p  +  p 

GD.  ..  =  C  ( | v.  ..i+c.  ..)  ':w - liii — 


'I'M  ~  ^2 M (pi;j+i  +  2'pi,7 +  pif j.p 


where 


Cr  =  0.02  £2  +  0.09  £2 

52  X(i,jj)  r(i,jj) 


c  =  0.10  n2  +  0.02  n2 

n2  r(i,jj)  x(i,jj) 


where  jj=j  and  j+1  for  the  predictor  and  corrector  steps, 
respectively  and  c^  ^  is  the  local  speed  of  sound. 

It  can  be  shown  that  this  damping  will  give  rise  to  fourth 
order  terms  in  the  form: 


C  AtA£2  ^  1  u !  +c 

^ ^  3£  [  4P 

3  2P 
3?2 

3ul 

+  C  AtAn  J- 
T)1  3r) 

|  V  I  +c 
4P 

3  2P 

3n2 

3U1 

3nJ 

in  the  ^-direction, 

and 

cn/“"3  A 

3  2P 

3^7 

au  1 

3n  _ 

3  9 

+  Cr  AtA£  ^ 

C2  as 

1  |n|+c 

L  4P 

3  2P 

352 

3u" 

3£_ 

in  the  n-direction . 

In  addition  to  this  pressure  damping,  the  normal  stress 

9 

damping  concept  as  used  by  McRea  and  others  was  used  here.  This 
damping,  represented  by  the  factor  3  in  Equation  (2.7),  affects 
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directly  the  normal  stresses  without  affecting  the  shear  stresses. 
The  value  of  3  =  -6  was  used.  It  is  not  always  possible  or 
recommended  to  use  very  large  values,  at  least  for  the  explicit 
schemes,  since  it  directly  decreases  the  allowable  time  steps, 
through  Afc,  as  can  be  noticed  from  Equation  (3.3).  It  is  preferable 
to  remove  3  when  convergence  is  almost  achieved  to  eliminate 
additional  changes  to  the  original  Navier-Stokes  equations. 

4 .  IMPLEMENTATION  OF  THE  BOUNDARY  CONDITIONS 

It  is  noted  here  that  the  steady  state  boundary  conditions 
for  the  time  dependent  equations  are  imposed  for  both  parts 
(Equation  (3.1))  of  the  original  governing  equations. 

First-order  accurate,  two-point  extrapolation  was  used  to 

compute  the  value  of  the  variables  at  the  wall,  far-field  boundary 

and  the  axis  of  symmetry  for  the  derivative-type  conditions. 

Although  arguments  may  arise  concerning  the  need  for  use  of  three- 

point,  second-order  accurate  forms  to  keep  the  numerical  scheme 
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totally  second-order  accurate,  it  was  experienced  that  this 
leads  to  instability  in  some  cases.  Further,  the  use  of  a  two- 
point  formula  is  more  consistent  with  MacCormack's  basic  two- 
point,  two-step  scheme. 

In  addition,  it  was  found  essential  during  the  early  time- 
evolvement  of  a  solution  to  fix  the  density  at  the  point  1=27 
(Figure  2.5b).  The  reason  for  that  is  that  the  pressure  at  point 
1=27  is  relatively  low  in  the  low  pressure  pocket  in  the  nozzle 
lip  region,  relative  to  the  value  at  point  1=28  (in  the  jet) . 
Therefore,  instead  of  fixing  the  pressure  at  this  point  (1=27)  , 
the  density  was  instead  fixed  to  influence  the  continuity  relation 
and  the  energy  equations,  thus  influencing  the  solution  more 
significantly.  The  average  value  for  points  1=26  and  28  was 
assigned  to  the  point  1=27.  This  procedure  was  maintained  until 
the  solution  stabilized  and  then  it  was  lifted  without  difficulty. 
However,  obtaining  a  solution  was  not  possible  without  such 
procedure,  due  to  blow  up  in  the  solution  due  to  negative  pressure 
values  occurring  in  this  low  pressure  region  during  early  stages  of 
convergence . 
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5. 


INITIAL  CONDITIONS 


Since  signals  can  only  affect  downstream  flow  computations  in 
supersonic  flow  cases  (along  lines  of  characterisitcs  and  within 
the  zone  of  influence) ,  it  was  anticipated  that  the  incoming  flow 
profiles  at  the  starting  station  will  strongly  affect  the  solution 
in  the  whole  domain.  Therefore,  the  initial  conditions  were 
based  on  that  profile. 

These  initial  values  are  not  the  most  adequate  near  the  jet 
exit  station,  but  the  boundary  conditions  at  this  location  are 
expected  at  the  supersonic  speed  of  the  jet  to  influence  the 
solution  in  this  location  immediately  in  time. 

It  was  noticed  that  starting  with  a  solution  that  allows  the 
flow  specified  by  the  incoming  profiles  to  expand  rather  than  to 
recompress  helps  in  obtaining  quicker  convergent  solution  and 
reduces  numerical  instability.  It  is  suggested  that  if  no 
appropriate  starting  solution  can  be  established,  then  starting 
with  low-level  values  everywhere  may  be  more  adequate  for  the 
supersonic  flow  cases. 
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SECTION  IV 

RESULTS  AND  DISCUSSION 


1.  DATA  OF  EXPERIMENTS  AND  THE  CASES  COMPUTED 

A  total  of  five  different  cases  were  computed.  Four  cases 
were  run  for  the  specific  nozzle  boattail  under  consideration  (the 
AGARD  10° -nozzle)  and  one  case  for  a  similar  AEDC  model,  with  cold 
jet  exhaust  to  study  the  effect  of  the  nozzle  boattail  geometry. 

The  experimental  results  of  the  base  case  with  the  cold 
exhaust  were  obtained  under  the  following  test  conditions  for  the 
external  flow: 

M  =  1.5,  T  =  580°R,  Re  =  2.5  x  106/ft 

oo  o  e  ' 

e  e 

with  the  jet  exhaust  being  under  the  conditions: 

TOj  =  54 0 °R,  RT  (Throat  Reynolds  Number)  =  1.773  x  10^, 

nozzle  exit  diameter  =  3.982  inches. 

These  conditions  for  the  jet  flow  will  yield  a  nozzle 
pressure  ratio  NPR  of  7.09,  where  NPR  is  defined  as  Po^/P^  . 
However,  the  ratio  P./P  is  only  3.746. 

3  °°e 

The  following  conditions  can  then  be  computed  for  the  external 

flow: 


U  =  1470.4  ft/sec 
°°e 

P  =  354.1  lb/ft2 

e 

Also,  for  the  jet  flow  one 
Uj  =  1040.6  ft/sec 

P  =  2510  lb/ft2 
o  . 


T  =  4  00  °R 

°°e 

,  PQ  =  1300  lb/ft2 

0 

can  compute: 

T.  =  450  °R 
3 

P_.  =  1324.6  lb/ ft2 
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The  hot  exhaust  case  was  computed  with  TQj  =  2100 °R  and 
the  hot  boattail  wall  case  was  computed  for  T  =  1100°R.  The 
exit  conditions  for  the  hot  exhaust  can  then  be  computed  as: 

U.  =  2052.1  ft/sec  ,  T.  =  1749. 4°R 

3  3 

Pq  =  2510  lb/ft2  ,  Pj  =  1324.6  lb/ft2 

Finally,  the  corresponding  2-D  case  and  AEDC  nozzle  model  case 
were  computed  for  cold  exhaust  temperature  as  in  the  base  case. 


2.  COMPUTATION  DETAILS 

All  the  computations  were  run  on  the  CDC-6600  computer  with 
a  mesh  of  30  x  39  points  in  the  r)“£  plane  with  Arj  =  0.03448  and 
AC  =  0.0263. 


Computations  were  expected  to  last  until  the  accumulative 


v-4 


physical  time  reaches  at  least  ^t^  where  t  ^  =  L/U^  =  12.9  x  10 

second,  where  L  is  the  length  of  the  physical  computation  region,  here 

as  26.0  inches.  The  convergence  criterion  was  based  on  the  condition 

that  the  maximum  change  in  the  static  pressure  over  one  t^  should 

be  less  than  five  percent.  Only  the  first  case  with  the  AEDC  nozzle 

geometry,  with  initial  conditions  as  prescribed  earlier,  was  left 

to  exceed  8t  ,  and  the  c  at  the  boattail  surface  was  monitored, 
ch  p 

The  values  of  c^  at  different  times  are  shown  in  Figure  (4.1). 


The  largest  changes  observed  were  usually  in  the  static 

pressure  and  particularly  along  the  center  line  r  =  0.  For  the 

base  case  of  cold  exhaust  for  the  AGARD  nozzle,  the  maximum  change 

in  the  pressure  along  the  center  line  was  4.5  percent  along  the 

time  interval  t  =  4  t  ,  to  t  =  5.9  t  ,  , 

ch  ch 

For  other  cases,  computations  were  stopped  after  about 
t  -  4  t  by  making  use  of  the  computed ' case  of  cold  exhaust  as 
the  starting  initial  conditions,  thus  also  eliminating  most  of  the 
numerical  difficulties  and  achieving  faster  convergence  rates. 


The  computation  time  for  the  typical  base  case  was  0.022 
sec/grid  point  per  cycle  per  one  t^.  Cycle  here  denotes  three 
sweeps:  two  in  the  redirection  and  one  in  the  ^-direction. 
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This  amounts  to  a  large  computer  time  and  it  is  mainly  due  to  the 
very  small  allowable  time  steps  for  the  explicit  scheme  as  required 
by  the  stability  condition  of  Equation  (3.3).  For  the  base  case, 

_  g 

this  maximum  time  step  was  0.26  x  10  sec. 

The  lengthy  computations  were  made  with  a  restart  feature, 
allowing  reasonable  computation  intervals  and  continuation 
according  to  the  computer  resources  available.  However,  implicit 
or  hybrid  schemes  are  strongly  recommended  for  study  and  possible 
successful  applications  to  the  present  flow  case. 

3.  COMPARISON  WITH  THE  EXPERIMENT 

The  computer  code  was  first  validated  by  comparing  with  the 
experimental  results  of  Reference  (6) ,  where  the  surface  pressure 
coefficient  is  provided  for  the  cold  exhaust  case  at  the  present 
Mach  number  of  1.5  . 

Figure  (4)  shows  the  computed  pressure  coefficient  c  = 

2  P 

(P  -  P^  )/(l/2poo  ),  compared  with  the  experimental  measure- 

e  ©  e 

ments  for  the  upper  and  lower  locations.  The  recompression  seems 
to  be  weaker  than  that  measured,  and  the  agreement  does  not  seem 
to  be  excellent,  resulting  in  estimated  pressure  drag  higher  than 
the  experimental  value.  It  can  also  be  seen  that  the  smaller  cp 
computed  near  x  =  0  is  mainly  due  to  the  imposed  pressure 
profile. 

The  main  reason  for  the  present  deviation  is  that  the  flow 

-S. 

pattern  ceases  to  be  axisymmetric  during  experimentation  due  to 

the  effect  of  the  model  support.  The  flow  turns  into  a  truly 

three-dimensional  flow  with  strong  swirling  cross  flow,  while 

computation  models  true  axisymmetry.  The  present  results  however, 

are  better  than  any  known  for  supersonic  flow  at  this  relatively 

high  nozzle  pressure  ratio.  The  results  of  Reference  (2),  for 

example,  agree  well  with  the  corresponding  experiment  only  for 

subsonic  flow  up  to  M^  =0.9  and  NPR  up  to  2.95  but  quickly 

deteriorate  for  supersonic  flow  at  M,*  =  1.1  and  NPR  =  4.77. 

It  is  also  expected  to  deteriorate  more  at  the  present  conditions 

of  M  =1.5  and  NPR  of  7.09  . 
oo  e 
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Details  of  the  flow  field  for  this  base  case  can  be  seen  in 
Figures  9A,  10A,  11A  and  12A.  From  the  Mach  line  contours  of 
Figure  9A,  it  can  be  seen  that  the  subsonic  region  protrudes 
considerably  down  stream  of  the  nozzle.  No  reversed  flow  region 
was  detected  for  this  case,  as  was  also  confirmed  by  the  experiment 
as  indicated  in  Reference  (6) .  In  fact,  no  reversed  flow  regions 
were  reported10  for  the  10°  and  15°  nozzle  tested  with  NPR  up  to 
15.8  . 

Another  useful  result  is  the  static  temperature  distribution 
which  is  essential  for  the  infrared  signature  analysis  of  the  jet 
exhaust  plume.  The  static  temperature  contours  are  shown  in 
Figure  (10A)  where  the  jet  temperature  drops  due  to  the  expansion, 
then  increases  when  the  recompression  starts.  It  is  also  noticable 
the  effect  of  the  discontinuity  of  the  surface  temperature  boundary 
condition  where  T  =  580°R  on  the  boattail  surface,  while  the  jet 
exhaust  static  temperature  is  450 °R. 

The  flow  pattern  of  the  exhaust  plume  is  of  the  repeated 
weak  reflection  system.  The  pressure  and  the  u-velocity  along 

j 

the  centerline  of  the  jet  are  shown  in  Figure  (4.2)  where  the 
values  are  normalized  by  the  value  at  the  starting  station  at 
i  =  1.  The  flow  expands  to  a  low  pressure  at  B  and  the  reflected 
expansion  wave  at  the  plume  surface  reflects  as  compression  wave 
(Figure  (1.1A))  causing  the  pressure  to  increase  in  the  region  BC. 
Then  it  starts  to  drop  in  region  CD  due  to  the  expansion  suffered 
from  the  reflected  shock  that  is  reflected  again  at  the  plume 
boundary  as  an  expansion  wave.  It  is  interesting  to  notice  the 
change  in  the  velocity  component  along  the  axis  of  symmetry. 

4.  BOATTAIL  GEOMETRY  EFFECTS 

The  effect  of  the  boattail  geometry  on  the  pressure  drag  is 
of  practical  concern  for  jet  engine  designers.  Optimizing  the 
boattail  geometry  to  give  minimum  drag  is  usually  a  major  factor. 
Therefore,  another  nozzle  geometry  case  (the  AEDC  nozzle)  was 
computed.  The  pressure  coefficients  for  both  geometries  are  shown 
in  Figure  (5) ,  where  the  AGARD  nozzle  shows  smaller  c^ .  However, 
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Figure  lib.  Axial  Variation  in  u  and  P. 
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when  computing  the  pressure  drag  for  both  the  geometries,  they 
were  almost  equal,  thus  giving  the  AGARD  nozzle  the  advantage  due 
to  the  larger  wetted  surface.  The  "kink"  observed  at  x  -  2.2  inches 
for  the  AEDC  model  lacks  reasonable  explanation. 

The  flow  pattern,  in  general,  was  similar  to  the  AGARD  nozzle 
case  with  no  major  differences. 

5.  HOT  EXHAUST  EFFECTS 

The  effect  of  the  hot  exhaust  is  propagated  upstream  to  the 
boattail  surface  through  the  subsonic  region  of  the  boundary 
layer.  The  pressure  coefficient  for  this  case  was-  smaller  than 
the  corresponding  cold  exhaust  case  as  can  be  seen  in  the  middle 
curves  of  Figure  (6) .  Although  the  difference  seems  to  be  small 
and  the  curves  seem  to  be  slightly  shifted  upward,  the  computed' 
pressure  drag  was  considerably  smaller  for  the  hot  exhaust  case 
by  the  ratio  of  79%  in  comparison  with  the  cold  exhaust  case. 

This  is  quite  important  if  we  consider  the  exhaust  temperature 
ratio  as  T0^  (hot)  /TQw.  (cold)  =  3.89. 

The  mixing  of  the  hot  jet  with  the  outer  stream  is  presented 
by  Figures  (11A,  12A) .  The  static  temperature  profiles  for  both 
the  cold  and  hot  exhausts  normalized  by  the  corresponding  jet 
static  temperature,  are  shown  in  Figure  (11A) .  The  symbol  j 
represents  the  point  location  index  on  the  axis  of  symmetry.  At 
the  location  j  =30,  which  marks  the  location  of  the  far-field 
boundary  (at  a  distance  of  four  times  the  jet  exit  diameter ) ,  the 
static  temperature  =  0.66  Tj (hot)  and  also  =  2.90  T^  .  For  the 
cold  exhaust,  the  temperature  at  that  location  is  0.77  T j  and 
=  0.87  Tw  .  It  is  obvious,  and  expected,  that  the  hot  exhaust 
plume  needs  longer  distance  to  attain  free  stream  temperature 
value . 

The  mixing  in  the  axial  velocity  component  is  shown  in 

Figure  (12A) .  At  the  same  location  of  j  =  30,  the  velocity  is 

u  =  1.15  Uj  and  =  1.61  for  the  hot  jet.  The  corresponding 

values  are  u  =  1.47  U.  and  =  1.04  U  for  the  cold  jet.  The 

j  ooe 
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same  remark  about  longer  distance  along  the  axis  of  symmetry  to 
attain  the  external  free  stream  conditions  holds  for  the  axial 
velocity  component  for  the  hot  exhaust. 

It  is  then  appropriate  to  observe  the  changes  in  the  velocity, 
temperature  and  pressure  along  the  axis  of  symmetry  for  this 
particular  case  of  practical  importance.  This  is  given  by 
Figure  (7) ,  where  it  is  clear  that  the  high  static  pressure  of 
the  jet  drops  quickly  to  match  the  "surrounding”  low  pressure 
condition.  However,  it  seems  that  the  axial  velocity  and  the 
static  temperature  will  take  longer  distances  to  get  in  equilibrium 
with  the  corresponding  external  flow  conditions. 

Finally,  the  density  contours  for  the  hot  exhaust  are  shown 
in  Figure  (8) .  It  was  always  noticed  that  the  velocity  components 
and  the  static  temperature  had  smooth  profiles,  while  the  density 
and  consequently  the  static  pressure  (computed  from  the  equation 
of  state,  using  p  and  T)  had  local  irregularities  in  their 
contours  and  profiles.  This  might  be  explained  as  a  consequence 
of  having  either  excess  or  less  than  enough  pressure  damping  as 
described  in  Section  (-3.3).  Another  explanation  is  that  the  density 

r 

is  computed  from  the  first  order  continuity  equation,  with  one 
boundary  condition  to  enforce.  In  addition,  the  transformation 
derivatives  were  not  totally  smooth  due  to  the  specific  coordinates 
used  and  the  many  requirements  imposed  on  them.  These  irregularities 
in  the  derivatives  are  expected  to  be  reflected  in  the  obtained 
solution,  however  this  does  not  explain  the  smoothness  in  the 
velocities.  Figure  (8)  is  a  good  example  for  the  irregularities 
in  the  density  contours,  which  also  bears  the  additional  effects 
due  to  the  linear  interpolation  used  for  contour  plotting  in 
addition  to  the  effect  of  large  grid  mesh. 

6.  BO ATT AIL  SURFACE  TEMPERATURE  EFFECTS 

In  addition  to  the  low  temperature  case  of  T^  =  580°R,  a 
case  with  T^  =  1100°R  was  computed  and  the  pressure  coefficients 
for  both  cases  are  displayed  at  the  top  of  Figure  (6) .  The 
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pressure  coefficient  is  much  smaller  for  the  hot  wall  case  and 
the  drag  is  67%  of  the  corresponding  cold  wall  value.  The 
irregularities  in  c  curve  near  x  *  10  inches  is  partially  due  to 
the  effect  of  the  discontinuity  in  "wall"  temperature  between  wall 
and  the  exiting  jet. 

It  is  of  interest  to  see  the  effect  of  the  wall  temperature 
on  the  static  temperature  field.  This  can  be  studied  from  Figure 
(10B)  where  the  changes,  as  might  be  expected,  are  only  confined 
to  the  surface  region.  Otherwise  ,  the  same  pattern  exists  as 
compared  with  the  cold  wall  case  (the  base  case)  of  Figure  (10A) . 

The  velocity  flow  field  for  this  case  was  plotted  in  Figure 
(13A)  where  the  direction  of  the  arrows  represents  the  flow 
directions.  However,  the  stream  lines  cannot  be  traced  properly 
by  the  arrows.  Therefore,  the  corresponding  stream  lines  we. 
computed  by  direct  integration  and  contours  of  the  streamlines 
are  shown  in  Figure  (13B) .  No  separation  bubbles  were  detected 
as  stated  before.  Also,  a  mass  conservation  was  made  to 
compute  the  diameter  of  the  jet  tube.  For  the  jet  exiting  from  an 

j 

exit  radius  of  3.982  inches,  conservation  of  mass  showed  that 
this  radius  increases  to  6.510  inches,  i.e.  163  %  for  this 
case. 

7  •  THE  TWO-DIMENSIONAL  CASE 

A  two-dimensional  model  for  the  AGARD  10° -nozzle  axisymmetric 

\ 

case  was  computed  by  minor  changes  in  the  code,  mainly  by  setting 
jQ  =  0  instead  of  1  in  the  Equation  (2.8).  The  pressure  coefficient 
is  plotted  at  the  bottom  of  Figure  (6)  with  comparison  with  the 
axisymmetric  case.  The  pressure  coefficient  for  this  case  is 
larger  than  the  axisymmetric  case  before  the  recompression  at 
x  -  6  inches  and  smaller  after  the  recompression.  It  is  surprising 
that  the  recompression  occurs  at  almost  the  same  location  for  all 
of  the  cases  considered.  This  location  seems  to  be  more  dependent 
on  the  NPR  which  is  constant  for  the  present  cases.  There  is 
also  unexplainable  irregularities  near  x  -  10.0"  for  the  two- 
dimensional  case.  However,  it  is  expected  that  the  disturbance 
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due  to  the  two-dimensional  geometry  to  be  larger  in  general  and  to 
extend  far  down  stream  longer  than  in  the  axisymmetric  case.  This 
can  be  best  shown  by  Figure  (9)  for  the  Mach  line  contours.  The 
sonic  line,  for  example,  protrudes  down  stream  in  the  two-dimensional 
case  further  than  in  the  axisymmetric  case. 
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SECTION  V 

SUMMARY  AND  RECOMMENDATIONS 


1.  SUMMARY 

It  has  been  shown  that  complex  flow  cases  can  now  be  computed 
as  a  step  toward  the  final  goal  of  computing  the  flow  past  a  real 
three-dimensional  aircraft  model.  Finite  difference  solutions 
for  supersonic  flow  at  Mach  number  1.5  past  axisymmetric  and  two- 
dimensional  nozzles  with  real  jet  exhaust,  were  obtained  using 
the  complete  Navier-Stokes  equations.  The  viscous  turbulent 
mixing  between  the  two  flow  streams  was  handled  through  the  present 
approach,  without  any  assumption  or  approximation. 

The  explicit,  time  dependent,  second  order  accurate  MacCormack's 
scheme  was  used  with  the  emphasis  being  on  the  demonstration  of  , 
the  success  of  the  present  unified  approach,  rather  than  on 
emphasizing  high  accuracy  or  computational  efficiency. 

20 

Five  different  cases  were  computed  and  reported  ,  all  at 
nozzle  pressure  ratio  PQ  ./Pc»  =  7.09  to  study  the  effects  of 

boattail  geometry,  jet  exhaust  temperature,  boattail  surface 
temperature  and  the  differences  between  the  axisymmetric  and  two- 
dimensional  cases.  The  surface  pressure  coefficient  and  the  pressure 
drag  were  the  main  subject  of  interest,  more  than  the  profile 
drag  (skin  friction),  due  to  the  unavailability  of  experimental 
results  for  the  latter  and  also  to  relieve  the  larger  grid  point 
number  and  the  longer  computation  time  that  is  needed.  Comparison 
with  the  experiment  is  reasonable  if  the  three-dimensionality 
of  the  experimental  results  is  considered. 

The  hot  exhaust  case  yielded  a  smaller  pressure  drag,  79% 
of  the  corresponding  drag  for  the  cold  exhaust  case.-  The  hot 
boattail  surface  case  yielded  a  smaller  pressure  drag,  67%  of 
the  same  corresponding . cold  wall  case.  The  computation,  assuming 
axisymmetry  gave  larger  pressure  drag  than  experimentally 
predicted.  The  location  of  the  recompression  on  the  boattail 
surface  was  not  changed  in  all  five  cases  computed.  It  is  assumed 
that  this  location  should  be  a  strong  function  of  the  nozzle 
pressure  ratio,  which  was  constant  in  all  the  cases. 
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2. 


RECOMMENDATIONS  FOR  FURTHER  STUDIES 


For  better  accuracy,  more  grid  points  will  be  needed  for 

computation  if  the  turbulent  profile  drag  is  to  be  estimated  with 

reasonable  accuracy.  This  will  directly  increase  the  computation 

time  to  impractical  figures,  therefore,  implicit  or  hybrid 
21  22 

schemes  '  should  be  studied  and  used,  if  possible,  to  avoid 
the  tight  stability  restriction  on  explicit  schemes. 

The  normal  shock  mode  (Mach  disk  mode)  for  the  present  flow 

case  has  not  yet  been  computed  successfully  for  the  higher  nozzle 

pressure  ratios.  Studies  should  be  made  toward  establishing  the 

factors  that  affect  such  computations  (e.g.  initial  conditions, 

23 

type  of  pressure  boundary  conditions  ,  mesh  spacing  near  and  across 

24 

the  disk,  different  approaches  or  numerical  methods  ,...).  Use 
of  the  strong  conservation  form  of  the  governing  equations 
(Equation  2.9)  rather  than  the  weak  form  (Equation  2.8)  might  be 
tested  to  access  the  effects  on  accuracy  and  possible  improvement 
in  shock  capturing  ability. 

The  effects  of  the  turbulence  eddy  viscosity  model  used  on 
the  obtained  solutions  have  not  yet  been  established.  Studies 
for  better  models  are  still  needed.  Also,  better  and  more 
representative  boundary  conditions  (profiles)  at  the  starting 
station  of  computation,  would  be  of  direct  help  for  improved 
results.  Such  information  can  be  obtained  from  experiments  if 
in-advance  coordination  between  the  experimental  and  computational 
efforts  are  possible.  In  addition,  the  effect  of  the  location  of 
the  far-field  conditions  on  the  solution  is  still  needed  to  be 
studied  and  evaluated. 
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DOMAIN  OF  PRESENT  INTEREST 


Figure  12 .  Domain  of  Present  Interest 


<:z 


IN  COMPUTATIONAL  P 
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Figure  14 .  Boundary  Conditions  in  the 
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sical  and  Computational  Domain. 


Figure  15.  Comparison  with  Experiment. 


Figure  16. 


Influence  of  Boattail  Geometry  on  Pressure  Coefficient 
Distribution . 
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Figure  18.  Axial  Variation  in  u,  p  &  T. 
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Figure  19.  Density  Contours. 


Figure  21.  Static  Temperature  Contours. 
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Figure  22.  Static  Temperature  Profiles  in  the  Mixing  Region 


Figure  23.  Axial-Velocity  Profiles  in  the  Mixing  Region 


Figure  24.  The  Flow  Field  and  Stream  Function. 
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